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We propose a general method for simplifying master equations by eliminating from the description 
rapidly evolving states. The physical recipe we impose is the suppression of these states and a 
renormalization of the rates of all the surviving states. In some cases, this decimation procedure 
can be analytically carried out and is consistent with other analytical approaches, such as in the 
problem of the random walk in a double-well potential. We discuss the application of our method to 
nontrivial examples: diffusion in a lattice with defects and a model of an enzymatic reaction outside 
the steady state regime. 



I. INTRODUCTION 

Many problems of relevance in physics, chemistry and 
biology are appropriately described as systems of in- 
teracting, discrete entities which evolve according to 
stochastic rules. These entities may be not all equal, 
for example they may be molecules of different chem- 
ical species. When a) there is spatial homogeneity and 
other continuous degrees of freedom are irrelevant for the 
description, and b) the Markovian property holds, such 
systems are well described by master equations, 

= {WranPm ~ WnmPn) ■ (1) 

m 

Here n denotes the state of the system, P„ {t) is the prob- 
ability of being in state n at time t, and Wmn is the 
transition rate from state m to state n. For instance in 
multi-species chemical systems, 71 is a vector (ni, . . . , Us), 
the component nj being the number of molecules of the 
j-th specie. 

The description of chemical and biological processes 
in terms of master equations is usually rather accurate. 
However, it can have rather severe problems, even in the 
numerical treatment, if the number of states is large and 
overall if many degrees of freedom are involved. When 
the number of particles is large enough, a possibility is to 
describe the dynamics in terms of concentrations: there 
are several methods to move from a master equation to 
a partial differential equation, like the Van Kampen sys- 
tem size expansion^. The most common result of such a 
procedure is a Fokker-Plank equation. 

On the other hand there are systems which are funda- 
mentally discrete and, therefore, the continuous approx- 
imation may give answers which are quite different from 
that of the discrete case, like enzymatic reactions which 
are often studied in the limit where one has many sub- 
strate molecules but very few enzymes^. Other notable 
examples are genetic systems working with low molecules 
copy numbers^i^ and ecological systems close to the ex- 
tinct absorbing stated. 

A common features of many of these problems which 
can be used to simplify the description is the presence 
of many relevant timescales. Sometimes the timescale 



of the interactions is fast, but due to the many degrees 
of freedom the timescale of the global dynamics is much 
slower. This is the case of protein folding: while the 
elementary time scale of vibration of covalent bonds is 
~ 10~^^s, the folding time for a protein may be of order 
of seconds. Another problematic situation is when the 
degrees of freedom are not so many, but the interactions 
are of diverse nature, resulting in entries in the transition 
matrix Wmn of very different orders of magnitude. An 
example of the latter case comes from many cell regula- 
tory systems: proteins in a cell can interact chemically, 
again on molecular timescales, and via transcription reg- 
ulation, on timescales of seconds or even minutes. In all 
these situations one says that the system has a multiscale 
character—. 

The necessity of treating the "slow dynamics" in terms 
of effective equations is both practical (even modern su- 
percomputers are not able to simulate all the relevant 
scales involved in certain difficult problems) and concep- 
tual: effective equations are able to catch some general 
features and to evidence key controls and basic ingredi- 
ents which can remain hidden in the detailed description. 
The study of multiscale problems has a long history in 
science: perhaps the first example is the study, due to 
Newton, of the precession of the equinoxes, which was 
basically a special version of the averaging method in 
mechanics^. In fluid dynamics an example of the multi- 
scale procedure is the derivation of an effective Fick equa- 
tion for the large scale and long time behavior starting 
from the transport equation^: the diffusion tensor de- 
pends, often in a non intuitive way, on the velocity field. 
Finally, in quantum mechanics, since the nuclei are much 
heavier than the electrons, one can simplify the treatment 
"splitting" the electronic and nuclear degrees of freedom, 
like in the Born-Oppenheimer approximation^, or in its 
modern generalization due to Car and Parrinello^. 

Multiscale systems described by master equations can 
be difficult to deal with. Gillespie's stochastic simula- 
tion algorithmii, which is an exact and well established 
method for numerical simulations of master equations, 
is not very efficient to simulate multiscale systems. The 
reason is that Gillespie algorithm treats fast and slow 
dynamics on equal footing, while if one is interested in 
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the slow dynamics it is usually not necessary (and com- 
putationally demanding) to exactly integrate the fast 
dynamics. In recent years, several authors introduced 
modifications and approximations of the Gillespie algo- 
rithm that allow to improve its efficiency when applied 
to multiscale problems. Many approaches have been de- 
veloped like the quasi-steady state approximatio n^^'^'^ , 
fast variables elimination method a^^'^^ , finite state pro- 
jection technique s'^' or continuous approximations of 
fast variablesiSii^. 

The method developed in this paper follows a differ- 
ent idea. We want to study master equation containing 
states with fast and slow dynamics (evolving with differ- 
ent characteristic times). Our approach is to write down 
an effective master equation describing the evolution of 
the slow states only. The physical recipe we are going to 
impose is that, looking at the system on a slow timescale, 
the time spent on the fast states can be neglected. These 
states are thus eliminated from the description and this 
brings to new effective transitions among the slow states. 
Clearly this approximation gives better results when the 
separation of timcscalcs between fast and slow states be- 
comes large. Our procedure, by reducing the number 
of states, may allow for a large gain in simulation time. 
However, our main goal is to obtain in a systematic way 
a simple description of a system obeying a master equa- 
tion on a slow timescale. Indeed, we will show that in 
some cases our method does not correspond only to a 
numerical recipe, but allows for analytical predictions. 

The plan of the paper is the following. In Sectionllllwe 
introduce our method. Section |llT] is devoted to the one 
dimensional examples: a random walk in a double- well 
potential and a random walk in a lattice with defects. 
Section IIVI wc apply our method to a common model 
of enzymatic reaction and show how it can predict its 
behavior far from the steady state regime. Conclusions 
and perspectives are in section [V] 
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FIG. 1: Sdiematic example: a particle confined in a four-well 
potential. 

If one is interested to properties on times much longer 
than Ti , it is quite natural to devise a model of the sys- 
tem with only two states, say A which includes 1 and 
2 and B which includes 3 and 4. This simple example 
suggests how the study of a system on a slower timescale 
may lead to a reduction in the number of states one has 
to consider. 

Consider now a generic master equation 

= {W^nPm - WnM ■ (2) 

The standard way to simulate the above equation 
would be of course the Gillespie algorithroii, which is 
exact and docs not imply a choice of a timestep. On the 
other hand, wc would like to select a (slow) timescale, so 
let us discuss what happens if we integrate naively the 
equation, for example discretizing the time via the Euler 
algorithm: 



II. THE METHOD 

Let us introduce the intuitive idea of the method with 
a schematic example. We consider the motion of a Brow- 
nian particle in a potential V{x) having 4 minima, like 
the one sketched in Fig. [T] 

When the temperature T is not too large the probabil- 
ity distribution function of the particle is concentrated 
only around the the minima xi,X2,X3 and X4. There- 
fore, instead of the complete Fokker-Planck description, 
it is sensible to study the system in terms of a Master 
equation with 4 states: we say that the system is in the 
state i = 1, 2, 3, 4 when x{t) is close to Xi. Moreover, we 
consider the case in which the barriers between states 1, 
2 and 3, 4 have the same height AVi, being smaller than 
the height AV2 of the barrier between states 2 and 3. In 
this case, one has two very different typical times: the 
transition time between states 1, 2 and 3, 4 ri ~ gAVi/T 
and the transition time between 2 and 3, T2 ~ e^^^^"^. 



m m 

(3) 

In this case, the application of the Euler algorithm 
corresponds to approximate the master equation with a 
Markov chain, defined by the Markov transition matrix 

4At ^ / At Wnm m ^ U , , 

^mn I I _ AtW°"* m = n ^ ' 

where we introduced for convenience of notation the total 
out-rate of state n. 

When integrating Eq.([21), At should be chosen to be 
small compared to the timescales of the system dynam- 
ics. However, notice that Eq.([21) and ([3]) share the same 
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stationary condition independently on At. This means 
that they have the same stationary state for any At, even 
though they can have in principle quite different dynam- 
ics when At increases. A natural question is what hap- 
pens when At becomes large. The problem one encoun- 
ters is that the diagonal element of the matrix ([4]) can 
become negative. This starts happening when 

At>max{{WrY'] ■ (6) 

On the other hand, when At < max„{(VF°"*)~i}, the 
Markov chain and the original master equation are closely 
related. Not only they share the same stationary state, 
but one can easily write a relation between their eigen- 
values, which are in a one to one correspondence. We 
can write the Markov matrix A^^*^ : 

^(At) ^^^t (7) 



A 




FIG. 2: Schematic example of the transition rates before and 
after the decimation in a master equation with 4 states. Here, 
the state C on the left is a fast states and is eliminated from 
the description. On the right we represent the renormalized 
transition in terms of the original rates. The new transition 
rate from B to D contains the original one plus an additional 
contribution coming from the elimination of C. The transition 
rate from A to D, equal to zero in the original graph, contains 
only the effect of the decimation. 



where / is the identity and W is the matrix having as 
non-diagonals elements the transition rates and on the 
diagonal minus the total out-rates. This implies that, 
calling Xa and Xw the eigenvalues of the matrices A and 
W, one has 

A^"^*^ = 1 + AtXw (8) 

Another way of seeing it is that A'^'^*^ is the first-order 
expansion of the evolution operator of the master equa- 
tion cxp{W'' At) . This implies the error one makes on the 
eigenvalues is order O(A^At^) and should be negligible 
for the slower modes when At is not too large. 



A. Decimation procedure 

Eq. ^ suggests a strategy to identify "fast states" 
on a given timescale. The idea is that, every time a 
generic state n is reached, the average time spent in it is 
As a consequence, if we choose the parameter 
At representing the smaller timescale we aim to describe, 
then all states n having W°^* > At~^ should not enter 
into the description. In order to eliminate them, the 
physical recipe we impose is simply that the time spent 
on these states is zero. In this way, the states disappear 
from the dynamics and the transition rates from a generic 
state k to a generic state j are modified according to: 

K = W,,+W,„^. (9) 

This procedure corresponds to adding to the rate the 
process of k going to n with the proper rate and then 
instantaneously going to j with the proper probability. 
A graphical example of the application of the method is 
shown in Fig. 



B. Commutativity and adiabatic approximations 

A natural question at this point is whether the or- 
der of elimination of the "fast states" matters for the 
resulting dynamics . In presence of two fast "linked" 
states, say n and to, such that At > max {(iy°"*)~^}. 
At > max{(V(C*)-i} and W^^ > or W„„ > 0, the 
result of decimation could be different, in principle, de- 
pending on the order of elimination of the "fast states" , 
i.e. before n and them m or viceversa. 

To show that our procedure commutes in general and 
clarify the connection with adiabatic approximations, let 
us rearrange the vector {Pi{t), P2{t) . . .) into two vectors 
{ipi{t),'ijj2{t)) where ^pl is a vector containing the proba- 
bilities of the slow states and ijj2 contains the probabilities 
of the fast states. We rewrite the master equation as: 

d_( Bn Bi2 \ ( ^i\^(Biii^i+Bi2^2\ 

dt\i^2 ) \ B21 B22 ) ) \B2li'l + B22^p2 J ■ 

(10) 

The adiabatic approximation corresponds to set 
dip2/dt ~ 0, that is 

B21^1 + S22V'2 = 0. (11) 

We can safely assume det 7^ since otherwise the 
probabilities of the slow states would be zero at equilib- 
rium. So we can solve the above equation for 'ip2 and 
substitute it into the equation for ipi: 

^7/-! = {Bu - Br2B-^'B2i)^i. (12) 
Notice that Eq. (fT^ preserves normalization since 

I ^[P„(t) e Mt)] - I E ^» W = (13) 

71 n 

in other words, within the adiabatic approximation there 
is no probability flow between fast and slow states. 
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Eq. ([12]) may thus be considered a bona fide reduced mas- 
ter equation for the slow degrees of freedom, with the 
transition rates being rcnormalizcd due to the effect of 
the fast states. Notice that, at variance with fast vari- 
ables elimination methods^^, the probabilities tpi are not 
marginalized probabilities, meaning that we didn't aver- 
age over fast states. The probabilities of fast states can 
be eventually reconstructed as a function of time after 
solving the equations for the slow states. 

To clarify the relationship between Eq. and our 

method, notice that condition (jlip is a linear set of equa- 
tions for the fast degrees of freedom. Let us consider 
the solution obtained by the substitution method: we 
start eliminating a particular fast state j by computing 
its probability from the j-th equation: 



P. 
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(14) 



It should be clear at this point that, by substituting 
the above expression in the evolution equation of a fast 
state, the original rule of Eq. ^ is retrieved; the same 
procedure can be iterated to solve the probabilities of all 
the fast states and substitute it into the remaining equa- 
tions. The conclusion is that our method is equivalent to 
solve the condition pT]) using the substitution method. 
But since the solution of that condition is unique, the 
resulting master equation will be Eq. p^ . independently 
on the order of elimination of the fast variables. 



III. ONE DIMENSIONAL EXAMPLES 

Now we apply our decimation procedure to some one 
dimensional random processes. The first example is the 
well known problem of the double well; in this case the 
decimation procedure can be carried out analytically and 
predict the correct transition rate between the two min- 
ima. The second example is a one dimensional random 
walk with defect. 



A. Double well potential 



FIG. 3: Particle in a discrete double well potential. A and B 
are the states corresponding to the two minima of the poten- 
tial. 



the probability of being in state n is proportional to 
exp(— /314,). We do not restrict oursclf to a particular 
choice of the potential, we just assume that the potential 
has two minima which we denote with n = A and n = B, 
as sketched in Fig. ([3]). As usual, we are free to choose 
a value of At and eliminate all the states whose inverse 
out-rate is smaller than Ai: 



At > l/W" 



(16) 



By increasing At, the two minima are the last two sur- 
viving states since in a minumum both the exponentials 
in Eq. (jl6p have negative arguments. This means that 
with a proper choice of the time scale we can end up 
with a two state Markov chain and calculate the transi- 
tion rate between the two minima. 

By ailing N the number of states between A and 
the transition rate writes: 
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(17) 



Let us consider the problem of a random walk in a 
double-well potential. This example is paradigmatic both 
in physics and in kinetic chemistry, in the latter case the 
one-dimcnsional axis represents the reaction coordinate 
and one is interested the rate of jumping from one mini- 
mum to the another, i.e. of the reaction to occur. Again, 
we assume that this axis is subdivided in discrete states 
{n}. We can introduce the potential Vn, which deter- 
mines the following transition probabilities: 



^j±i = exp 



f3 



(15) 



being f3 the usual Boltzmann factor. Notice that 
the above transition rates ensure that at equilibrium 



In the denominator we indicate with the notation 
the total out-rate of state i to remember that this out- 
rate changes when a neighboring state is eliminated. This 
means that the product of out-rates has to be evaluated 
by eliminating the states one after the other. To evaluate 
the expression ()17p . we make use of the following equality, 
which is demonstrated in the Appendix: 

i+N N / i+] i+N 

n = E n ^fc-^-i n ^^-^-+1 

h=i+l j=0 yk=i+l k=i+j + l 

(18) 

with the convention that products of less than one terms 
arc always equal to one, Hfc^i Wk = l. Substituting H 



and (1181) into (fTTI) one obtains 



eM^{VA'VB)] 



Ef=o eMi(yA+,~VA + VA+j+i-VB)] 
1 1 



= (19) 



N N 

E exp[|(KiH.,+F^+,+i)] E exp(/?V:4+,) 

where in the last step we assumed that the potential does 
not change much between adjacent states. The last ex- 
pression is basically the well known result for the transi- 
tion rate in a double well potential in a continuous system 
described by a Fokker-Planck^S. 
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B. Random walk with defects 

Considered now a random walker in a one dimensional 
geometry with periodic boundary conditions; the number 
of its possible states is fixed to TV = 10'^. The transition 
rate between two adjacent states is Wn.n-\ = ^n,n+\ = 
1 except for a fraction p = 0.2 of the states (the "de- 
fects"), having Wn,n-i = Wn,n+i = 20. The defects are 
placed at random at the beginning and kept the same 
in all the simulations (quenched disorder). In solid state 
physics the master equation of this system corresponds 
to the Schrodinger equation (at imaginary time and with 
disorder) in the limit of tight binding approximation^!. 

Clearly one has 1^°"* = 2 for the normal states and 
y^out ^ j-Qj, ^j^g defects. This means that the defects 
are the fast states and we will test what happens when 
eliminating them from the description. For comparison, 
we will see also the differences with the ideal case without 
defects, that is Wn,n-i = Wn,n+i = 1 for all the states. 

It is easy to analytically calculate the eigenvalues in 
the case without defects^i: 



A, 



1 



i = 0,±1,±2. 



N 



1 



(20) 

Notice that all the eigenfunction are non-localized, 
which in solid state physics corresponds to a conductive 
system. On the contrary, in the case with defects (for 
any positive value of p) the eigenfunctions are localized, 
i.e. one has a transition from a metal to an insulator—. 

In Fig. (HI we show the eigenvalues of the master equa- 
tion in the three cases: random walk without defects, 
with defects and after decimation. Being the master 
equation linear, the eigenvalues spectrum contain all the 
informations about the dynamics. The spectrum shows 
an interesting property: the most negative eigenvalues 
are much larger in modulus and correspond to fast de- 
caying eigenfunctions concentrated on the defects. It is 
clear from the figure that the effect of our algorithm is to 
eliminate these eigenvalues (and the corresponding eigen- 
functions). 

It is clear in the figure that the 100 eigenvalues smaller 
in modulus, corresponding to the slower dynamics of the 



FIG. 4: Eigenvalues, listed in decreasing order of a random 
walk with defect. The number of state is A'' = 10"^. Black 
circles: no defects, the jump rate is Wi.i+i — Wi.i=i = 1. 
Red squares: with probability p = .2 the sites have defects 
and their jump probability is VPTji+i ~ W^m=i = 20. Green 
diamonds: random walk with defects after applying the dec- 
imation scheme to all the fast states. All the eigenvalues are 
real. The difference between the figure and the inset is just 
the axis scale. Notice in the main figure that the first 100 
eigenvalues of the decimated problem follow very closely the 
case with the defects. In the inset the eigenvalues correspond- 
ing to the fast states are evident; notice that their number is 
the same as the average number of defects p * N. 



system, are very similar in the model with defects and the 
decimated one. A further check comes from the correla- 
tion functions. We performed the simulations starting 
from a slow state i = 500 and plotted in Fig. [5] the prob- 
ability of being in the state i as a function of time in the 
three cases we considered. After an initial time, the cor- 
relation functions for the decimated model and the model 
with defects are very close to each other. 



IV. ENZYMATIC REACTIONS 

In this section we apply our method to a model of en- 
zymatic reactions. Enzymatic reactions are widely stud- 
ied in kinetic chemistry and many methods have been 
derived to predict the steady-state kinetics in cases in 
which the reaction is close to an equilibrium, or at least 
in a non equilibrium steady state^^. To exemplify the rel- 
evance of our method in this case, we will consider in the 
following the simplest model of an enzymatic reaction: 



S ■ 



ES 



E 



(21) 



where as usual E denotes the enzyme, S the substrate 
and P the reaction product. This correspond to the fol- 
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the out-ratcs suggests two limiting cases in which there 
is separation of scales between fast and slow states. The 
first case is kiNs <C + k2). In this case, states hav- 
ing Ne > have a much greater out-rate than those with 
Ne = and can be coarse-grained. Physically it corre- 
sponds to the situation in which the complex is efficiently 
processed, so that its concentration is always very close 
to zero. In this case the application of our method yields: 



Krn. 



(25) 



Notice that we wrote the rate in term of N'g (the total 
number of substrate molecules) and not Ns since with the 
coarse graining we grouped together states with bounded 
and free substrates as long as the total number is the 
same. Notice also that this rate correctly corresponds 
to the Michaelis Menten formula, eq. ([M]) taken in the 
limit Ns / Km ^ 1. The new information provided by our 
method is that in this limit the linear dependence on the 
number of substrates is valid also beyond the steady state 
approximation, that is when the enzymes concentration 
becomes large. 

The second limiting case is kiNs ^ {k-i +k2)- In this 
case, all states having Ne > and Ns > are fast. This 
means that, for a fixed N'S, the slow state is characterized 
by Ne = max(0, N'^-NJ). When N]^ < TVj we retrieve 
Ne) + PiNs + 1, Ne + 1) + (22) again the steady state result: 



FIG. 5: Correlation functions. The system is the same of 
Fig.Q prepared in the initial state i — 500 and the prob- 
ability of being in state i, averaged over 10^ realization, is 
plotted as a function of time. The three curves are: (black, 
continuous) the system without defects, (red, dashed) system 
with defects, and (green, dot-dashed) the decimated system. 



lowing master equation: 



-PiNs,NE) = ki{Ns - l)iNE - l)PiNs~l,NE-l) 



+k2iN'E~NE)P{Ns, Ne + 1) - P{Ns, Ar£;)M^°"*(7V£;, Ns) 

where we called Ns the number of free (non bindcd) 
substrate molecules, Ne the number of free enzymes 
molecules, and N^ the total number of enzyme 
molecules. We also introduced, as usual, the total out- 
rate of a general state: 



(26) 



W°''\Ne,Ns) = kiNENs + [k^i + k2){Nl 



Nf 



The above expression, again, is consistent with the 
Michaelis-Menteen equation, this time in the limit 
Ns/Km <C 1 when the reaction rate does not depend 
anymore on the substrate concentration. On the other 
hand, when iVj > A^J one has: 



(A:_ 



k2)Nl + NE[kiNs - (fci + k2)]. (23) 



W. 



koNi 



(27) 



The quantity of interest is usually the production rate, 
that is the velocity of the rightmost reaction as a func- 
tion of the concentration of enzymes and substrate. It is 
easy to estimate it in the quasi-equilibrium approxima- 
tion, that is when ^ k2, meaning that the leftmost 
part of the reaction can be considered at equilibrium. 
Another thing one can assume is steady state dynam- 
ics: the calculation is simple even if the complex is not 
at equilibrium but the concentration of the complex ES 
does not vary much with time, d[ES]/dt « 0. The steady 
state approximation is known to hold very well when the 
concentration of enzymes is much smaller than that of 
substrate. 

In this case the reaction velocity follows the Michaelis- 
Menten formula: 



k2Ej^S 



K 

where K,n = (fc-i + ^2)/^! 
method, we start from Eq. 



S 



(24) 



In order to apply our 
The expression for 



This last case is radically different with the steady state 
prediction: we have a linear dependence on the number 
of substrate molecules while Michaelis-Menteen formula 
predicts a rate which is independent on iVj. The phys- 
ical reason is that this is a case in which the complex 
never reaches a steady state: on a fast timescale, order 
fcf ^ the free substrate is completely converted into the 
complex ES. Then it is depleted on much slower times, 
and the depletion speed is limited by the lower between 
the enzyme and the substrate concentration. 

In Fig. © we plot simulations of the reaction (PT]) 
in four different situations. The top left figure cor- 
responds to the situations in which iVj <C A^J and 
kiNs ^ {k-i + fe) (numerical values of the parameters 
are in the figure caption), together with the prediction 
of Eq. , while on the top-right we show iVj <C N]^ 
and kiNs 3> (A:_i -I- ^2) with the prediction of Eq. ((26|) . 
These two are the cases in which steady state approaches 
are known to work and our approach simply corresponds 
to the limiting cases of Michaclis-Mcnten formula. On 





FIG. 6: Simulations of the master equation corresponding to 
the enzymatic reaction of Ea. (|2H) . In all simulations we start 
with Ns = 100 molecules of free substrate and let the system 
evolve; the rate as a function of the number of molecules is 
evaluated by averaging over 10* realizations of the process. 
The number of enzymes is N'£ = 1 in the top figures and 
iVU — 100 in the bottom figures. The reaction rates are (left) 
ki — 10~^, k-i = k2 = 1 and (right) fci = 1, k-i = k2 — 
W~^. The lines in the black and white versions of the right 
figures are barely visible since the points fall very close to 
them. 



the bottom figures we show the same cases, but with 
N"^ < Ng . Notice that the left one correspond to the 
top/left one, rescaled with the higher number of enzymes; 
in particular, it is still described by Eq. (|25|1 . On the 
other hand, the bottom-right is even qualitatively differ- 
ent from the corresponding top figure. Eq. (j27p correctly 
describes the behavior, while Michachs Mcnten formula 
would predict the rate to be independent of iVj. 

In order to study also the fluctuations of the process, 
we simulated both the full master equation and the sim- 
plified processes defined by Eq. ^ and (gT]). Then 
we compare in Fig. ([7]) the variance of the number of 
product molecules as a function of time. The result of 
the simulations is that the fluctuations are almost in- 
distinguishable in the reduced processes: the intuitive 
reason is that most of the fluctuations are related to the 
slow states and thus unaffected by the coarse graining 
procedure. 

We conclude this section with some remarks. The con- 
dition kiNs ^ (fc-i + fe) has a physical interpreta- 
tion: it corresponds to the situation in which the con- 
centrations are very low. Indeed ki is the only rate 
that depends on the volume since it is essentially de- 
termined by the search time. Conversely, the condition 
kiNs 3> (fc_i + k2) corresponds to large concentrations. 
This means that all the four cases we discuss are, in gen- 
eral, experimentally accessible. This kind of approach 
may be useful in the study of enzymatic reactions in- 
side the cells. Such reactions are often characterized by 
low number of molecules^, and probably the Michaclis- 
Menten picture is appropriate even with a low number of 



FIG. 7: Product variance as a function of time, each aver- 
aged over 10^ realizations. The four figures correspond to 
the same parameter choices of Fig. ([6]). Black continuous 
lines correspond to simulation of the full master equation, red 
dashed lines are simulations of the reduced processes defined 
by Eq.lHSJ, dHJ and ([^Tj). 



substrate molecules. On the other hand, one can study 
situations in which concentrations are very high, or cases 
in which the changes that can occur in the cell environ- 
ment, for example as a response to a rapidly varying ex- 
ternal signal^!, can bring the reaction outside the steady 
state regime and make a description of this kind more 
appropriate. 



V. CONCLUSIONS 

In this paper we introduced a general method for the 
decimation of "fast modes" in systems evolving accord- 
ing to a Master equation via a coarse graining procedure. 
Our method is general, and somehow, in the spirit of the 
rcnormalization group (RG) approach. Indeed, similar 
approaches have been proposed in statistical mechan- 
ics for the study of disordered systems^^. At variance 
with the RG, we do not aim at reaching a fixed point 
by repeating the decimation procedure. In general, the 
method is applied only once: a sensible value of the pa- 
rameter At, which selects the coarse graining level, may 
be easily chosen by looking at the magnitude of the total 
out rates W^°"*. 

We show that the procedure is commutative and brings 
to consistent results for a general master equation, inde- 
pendently of the dimensionality. However, for the sake of 
simplicity, we discussed in details low dimensional exam- 
ples, where the method brings also to analytical predic- 
tions. In the discrete version of the double well potential, 
the decimation procedure is able to reproduce the well 
known result for the transition time. A numerical analy- 
sis shows that the method gives a very good approxima- 
tion of the original system in the case of a random walk 
with defects. In the case of an enzymatic reaction, we 
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show that the method allows for predictions only in some 
well defined limiting cases. However, in these cases, the 
prediction arc more general than those obtained within 
a steady-state approximation. 

Let us conclude with a short comparison with other 
approaches to multiscale master equations. Our method 
is more similar in spirit to projection method a^^'^^ than 
to other quasi-steady states mcthod o^^i^^ . The reason is 
that the result of our procedure is still a master equation 
while the other method describe the fast dynamics in a 
different way (with a differential or Langevin equation). 
The main difference between our method and the finite 
state projection is that our procedure is "local" (we con- 
sider fast and slow states) while projection methods gen- 
erally consider eigenvalues and eigenvector of the tran- 
sition matrix. Projection methods are usually of more 
general applicability, however our procedure may allow a 
more transparent interpretation of the surviving states, 
as we show in the examples we considered. 

On the other hand, fast variables elimination 
methodsi^ii^ aim at writing an evolution equation for 
the probability of slow states summed over the probabil- 
ity of the fast ones. In this way one loses informations 
about the fast states and, if there is separation of scales, 
may write a closed equation for the slow ones. We show 
that our approach does not imply such coarse graining. 
In fact, no information about the fast states is lost in 
our case and their dynamics can be reconstructed after 
solving the problem involving the slow states only. 

There are also analogies between our method and some 
coarse graining procedure that have been proposed in the 
field of complex networks. In this field, a relevant prob- 
lem is the identification of communities2^ and a possible 
way to do it is to consider a diffusion process on the net- 
work and try to coarse grain the graph while keeping the 
long timescale properties of this dynamics^li^. The dif- 
ference is that in these cases the procedure allows for a 
spatial simplification of the links, while our decimation 
procedure in a system whose states are seen as elements 
of a graph brings a suppression of the fast states but 
without a relevant simplification of the surviving connec- 



tions. 



Appendix: proof of equality (|18p and commutativity 

in ID 



In this appendix we prove that, when decimating a 
cluster of N consecutive states, the product of their out 
rates can be written as: 

i+N N / i+j i+N 

n = E n ^k^k-i n ^^-.^+1 

h=i+l j=0 yk=i+l k=i+j + l 

(28) 

remembering the convention that the product of less than 
one object is equal to one, and it is independent of the 
order of decimation. 

i i+N+1 

•■■■O'O'O'O'O'O'O'OO'O'O'O'O'O'O'O'O'O* 



FIG. 8: Two states (i and i + N + 1, black) separates by a 
cluster of N states (white) to be decimated. 



We will show that the above formula holds when dec- 
imating states in consecutive order, say from state i + 1 
one after the other to state i + N (see Fig.®), remem- 
bering that the result does not depend on the order of 
decimation due to the commutative property. We will 
demonstrated it by induction: first of all, the above for- 
mula is obviously true for N = 1: 



W^,+i_,+2 + W,+i^^ (29) 

Now we show that if the equality holds for a given value 
of N, then it holds also for N + 1: 




i+N+1 N / i+j i+N \ 

II wr-" = E n ^k^k+. n ^Tn+i = 

h=i+l j=0 yk=i+l k=i+j + l J 

N / 1+] i+N \ / yri+N+l \ 

= E n w^-^^^ n ^^-^^^ w^.+.+M+^+2 + t:;-^ = 

j=o yfe=i+i fe=4-i-j-i-i / \ iifc=j+i ^k J 

N I i+j i+N+1 \ i+N+1 N+1 I i+j i+N+1 

= E n n ^k^k+i + n ^^^^-i = e n ^k^k+^ n ^^-^+1 

j=0 yA:=j+l k=i+j+l I k=i+l j=0 yk=i+l k=i+j+l , 

I 



which is exactly the same expression of eq. (pS)) for N 



+ 1 . This completes our proof. 
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